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The problem of packing a system of particles as densely as possible is foundational in the field 
of discrete geometry and is a powerful model in the material and biological sciences. As packing 
problems retreat from the reach of solution by analytic constructions, the importance of an efficient 
numerical method for conducting de novo (from-scratch) searches for dense packings becomes crucial. 
In this paper, we use the divide and concur framework to develop a general search method for 
the solution of periodic constraint problems, and we apply it to the discovery of dense periodic 
packings. An important feature of the method is the integration of the unit cell parameters with 
the other packing variables in the definition of the configuration space. The method we present led 
to improvements in the densest-known tetrahedron packing which are reported in Ref. Ij . Here, we 
use the method to reproduce the densest known lattice sphere packings and the best known lattice 
kissing arrangements in up to 14 and 11 dimensions respectively (the first such numerical evidence 
for their optimality in some of these dimensions). For non-spherical particles, we report a new dense 
packing of regular four-dimensional simplices with density (j) = 128/219 ~ 0.5845 and with a similar 
structure to the densest known tetrahedron packing. 

PACS numbers: 61.50.Ah, 45.70.-n, 02.70.-c 



I. INTRODUCTION 

The dense packing behavior of a general solid body 
(particle) in a Euclidean space is a problem of interest 
in mathematics, physics, and many other fields. A pack- 
ing is a collection of particles in the Euclidean space K'*, 
wherein no two particles overlap (i.e., the intersection of 
any two particles has an empty interior) and the pack- 
ing fraction or density is then the volume fraction of 
space covered by the particles. Of particular interest are 
packings of a given particle (wherein all particles are con- 
gruent), and the problem of interest is to determine the 
maximum possible density 0niax among all packings of 
a given particle. A packing that realizes this maximum 
can be thought of as the equilibrium state of the system 
of classical hard particles in the limit of infinite pressure 
or zero temperature. 

The general problem of packing congruent particles 
was posed as a part of the eighteenth of David Hilbert's 
famous Mathematische Probleme: 

How can one arrange most densely in space 
an infinite number of equal solids of a given 
form, e.g., spheres with given radii or regular 
tetrahedra with given edges (or in prescribed 
position), that is, how can one so fit them 
together that the ratio of the filled to the un- 
filled space may be as large as possible? [1] 

This part of the problem has been taken over the years as 



the resolution of the Kepler conjecture about the dens- 
est packing of spheres in three dimensions Q, and has 
therefore been considered resolved since the latter was 
proved by Hales Q- However, Hilbert's statement of 
the problem does not single out the sphere, and actu- 
ally mentions the regular tetrahedron as another parti- 
cle of interest. Recent work diverging from the focus on 
spherical particles has spotlighted ellipsoids regular 
and semi-regular polyhedra fw, 't'I (and the regular tetra- 
hedron in particular ^li, JTj ) , and superballs jl^]. Few 
bounds are known for the maximum packing fraction of 
general convex particles. Kuperberg and Kuperberg have 
shown that for any convex particle in two dimensions, 
(/>max > V3/2 w 0.86602 [H]. Torquato et al. used the 
known maximal packing density of spheres to derive an 
upper bound on the packing density of any solid, but 
this bound is trivial (i.e., 0max £ <t>^ , where cjP > 1) for 
many solids 6]. Ulam has conjectured that in three di- 
mensions, the sphere achieves the lowest maximum pack- 
ing fraction, 0i„ax = 7r/-\/T8 ~ 0.74048, among all convex 



particles [14 1 
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In the quest for dense packings of various particles, 
analytic and numerical investigations have both played 
important roles. The former have been very successful 
in the study of the dense packing of spheres, where an- 
alytic constructions based on groups, codes, and lami- 
nated lattices have produced the densest-known sphere 
packings and lattice sphere packings in many dimensions 
[15| . However, the analytic approach to the construction 
of dense packings relies on the imagination of the con- 
structor, and for a variety of other problems the densest 
packings have evaded the creativity of analytic investiga- 
tors and were only uncovered in computational investiga- 



2 



tions. While complete (i.e., exhaustive) algorithms exist 
for some problems (such as the algorithm in Ref. [l6j . 
which gave new best known results for the lattice cover- 
ing and covering-packing problems in some dimensions) , 
they do not exist or have runtimes that are too long for 
other problems. In those cases, incomplete search algo- 
rithms become necessary. 

One example of a dense packing that has only been 
uncovered by a de novo numerical search is the cur- 
rently densest-known packing of tetrahedra, whose struc- 
ture was first hinted at by a numerical search using the 
method described in this paper [ij. The structure was 
later optimized by Torquato and Jiao [lOj and by Chen 



et al. Results of subsequent Monte Carlo simu- 



lations have reproduced this structure and suggest it is 
the densest packing of regular tetrahedra at least with a 
small number (< 16) of tetrahedra in the unit cell (lol . [Tl| . 
Another de novo search with Monte Carlo dynamics has 
uncovered a packing based on a quasicrystal approximant 
reminiscent of the Frank-Kasper cr-phase with a slightly 
lower density foj . As these two structures were overlooked 
by previous analytical investigations [1, it is quite 
likely that without the results of de novo searches, they 
would have remained unimagined and undiscovered. 

In the best case, such searches would produce the op- 
timal packing possible subject to the built-in restrictions 
(such as number of particles in the unit cell or unit cell 
shape). However, in problems exhibiting a large degree 
of frustration, the presence of many local optima that are 
separated from each other by high barriers complicates 
the task of finding the optimal packing. The tendency 
of simulations to get stuck in the local optima of such a 
rugged optimization landscape, especially when these lo- 
cal optima proliferate as more particles are simulated, has 
been held responsible for suboptimal results in searches 
[1, 0|- technique which has been observed to re- 

lieve dynamical stagnation in Monte Carlo simulations 
at high pressures has been to allow slightly unphysical 
moves, such as allowing particles to temporarily overlap 
i- 

We propose a novel search method as an alternative to 
Monte Carlo simulations, with a number of features that 
directly address these observations. The method is based 
on the dynamics of the difference map, a constraint- 
satisfaction iterative search algorithm, and on the di- 
vide and concur constraint framework (we abbreviate 
this combination D — C. where the minus sign stands 
for the difference map) |18l - [20| . It adapts the D — C 
approach to the case of periodic problems and we shall 
call it periodic divide and concur (PDC). The difference 
map is designed to avoid being trapped in local optima 
and has been demonstrated in multiple applications to 
find solutions of highly non-convex problems, including 
finite packing problems with large numbers of particles, 
from random starting configurations [T8l - [20j . The search 
proceeds through a non-physical configuration space, cut- 
ting through the conventional physical optimization land- 
scape. Still, it is to be expected that the exponential 



growth in the number of local optima in the configura- 
tion space, which the search will still have to traverse, will 
nevertheless lead to suboptimal results when many inde- 
pendent particles are included in the search. Therefore, 
as discussed below, it is crucial for the unit cell variables 
to be aggressively optimized so that the number of par- 
ticles to be simulated can be reduced. The incorporation 
of the unit cell variables directly into the basic dynamics 
of the search achieves this goal. 

Numerical searches are restricted to finite-dimensional 
configuration spaces, and therefore have been largely lim- 
ited to investigating periodic packings, packings which 
are preserved under translations by a lattice A. In a gen- 
eral periodic packing, the particles are partitioned into 
p orbits of the lattice A, and when p — 1 the packing is 
called a lattice packing. In physics, any periodic arrange- 
ment is usually referred to as a lattice and the special 
case of p = 1 is known as a Bravais lattice. In general, 
the maximum density need not be realizable by a periodic 
packing, but arbitrarily close densities are realizable with 
periodic packings of arbitrarily large p. Similarly, arbi- 
trarily accurate approximations of any packing can be 
obtained using a sufficiently large cubic or orthorhombic 
unit cell. However, due to the rapid increase in compu- 
tational complexity and the proliferation of local optima 
as the number of independent particles rises, it is of- 
ten preferable to include fewer particles but allow for a 
variable unit cell shape. We focus then on searching for 
packings with a small number of particles in the unit cell. 

To our knowledge, variable unit cells have only been 
introduced recently to searches for dense packings, for 
instance with the adaptive shrinking cell scheme in Refs. 

0] and with the use of Parrinello- Rahman dynamics in 
the space of lattices in Ref. 21]. The increased particle 
population associated with restricting unit cell variabil- 
ity can sometimes be tolerated in two and three dimen- 
sions, but in high dimensions the number of particles that 
must be simulated grows exponentially due to the curse 
of dimensionality and this approach becomes impracti- 
cal. The constraint-satisfaction formulation of the peri- 
odic packing problem used in PDC features a variable 
unit cell and naturally treats the positions of particles 
in the unit cell and the unit cell parameters on the same 
footing. This new approach allows us to successfully look 
for dense sphere packings in dimensions as high as 14, 
further than probed by any previously reported unbiased 
numerical exploration of periodic packings. 

Besides the density of a packing, another attribute of 
interest is the coordination number, that is, the num- 
ber of nearest neighbors of particles in the packing. In 
the case of spherical particles, this amounts to the num- 
ber of spheres in contact with a given sphere, known as 
the kissing number [15] . Searching for high-coordination 
number arrangements around a single sphere has been 
accomplished previously with the D — C method [l9[. 
Here we apply PDC to search for space-filling periodic 
arrangements of high coordination number, and particu- 
larly lattice arrangements. 
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An efBcient de novo numerical search method can pro- 
vide critical utility in the field of packing. In addition to 
the ability of a de novo search to provide confidence in a 
putative, but not proven, optimal result, a de novo search 
has often been responsible for surprising new results: two 
recent examples in which unexpected (as it turns out, 
quasiperiodic) packings were found as the results of de 
novo searches are in the problem of tetrahedron packing 
g| and in the ten-dimensional kissing number problem 
23 ]. It is with these motivations that we introduce the 
PDC method in this paper. In Section [ll] we introduce 
the D — C scheme by presenting a simple example which 
serves to motivate the constructions in the subsequent 
sections. In Section ITlll we formulate the problems tack- 
led in this paper — sphere packing, the lattice kissing 
number, and polytope packing — in terms of constraint 
satisfaction. In Section HVl we describe in detail aspects 
of our implementation of the PDC search, including ef- 
ficient computation of projections to the constraints of 
Section Hm In Section|V]we present some results of PDC 
for the problems discussed, including a newly discovered 
packing of regular four-dimensional simplices. In Section 
IVII we present concluding remarks. 



II. MOTIVATION 

A. The D -C scheme 

The key step in applying the D — C approach to pack- 
ing problems is to recast the problem as a problem of 
constraint satisfaction. Particularly, we must express it 
as the problem of finding a configuration in a Euclidean 
configuration space (17), which satisfies two constraints. 
We identify a constraint C with the subset C C of 
configurations satisfying the constraint. A projection of 
a configuration a; to a constraint C is the operation of 
finding a configuration x' € C that minimizes the dis- 
tance \\x — x'\\. Each of the two constraints {C,D C $7), 
must be simple enough that the operation of projecting 
an arbitrary configuration to it can be computed effi- 
ciently. The iterative map used in exploring the config- 
uration space takes advantage of the formulation of the 
problem in terms of two simple constraints, as outlined 
in section fllBI In this section we present the application 
of the D — C scheme to finite sphere packing problems, 
which has been developed and implemented in Ref . [l^ , 
as an introduction to the main ideas of the scheme. 

The defining constraint of packing problems is the con- 
straint that no particles in the packing overlap, which we 
call the exclusion constraint. As a simple illustration of 
this constraint, consider the exclusion of a pair of unit- 
radius disks in R^. In this case, the configuration space 
^' is parameterized by the positions of the centers of the 
two disks: 



The exclusion constraint is then 

i^cxci = {(xi,X2) e *: ||xi -X2II > 2} C (2) 

This constraint adheres to the simplicity criterion of hav- 
ing an efficient method to compute a projection to it. 
Specifically, the projection is given by 



TTifexcl [(X1,X2)] 



(x'i,xy if IIX1-X2II <2 
(xi,X2) otherwise. 



where, 



=Xl 



X2 =X2 - 



Xi - X2 



2||xi - X2II 
2- ||xi -X2I 
2||xi - X2II 



(xi 



X2j 



(Xl - X2) 



(3) 

(4) 
(5) 



as illustrated in Figure [T] 

A more complicated case arises when three or more 
disks are considered. In this case, the exclusion con- 
straint, 

K^^ci = {(xi,...x„) e «-: 

I |xi — Xj 1 1 > 2 for all 1 < z < j < 

is not a simple constraint according to the criterion 
above. Alternatively, we could replace iiTexci by many 
pairwise exclusion constraints 



{(xi,...x„) e ||x. -Xjll >2}. (6) 



The pairwise constraints are all individually simple. 
However, as noted above, we are limited to problems de- 
scribed by only two simple constraints. 

Divide and concur provides a general procedure for re- 
ducing the number of simple constraints to two, at the 
expense of enlarging the configuration space. This re- 
duction is achieved by parameterizing the configuration 
space with more variables than are necessary to fully 
specify a configuration. In the example at hand, the new 
configuration space is 

n {(xi^2,-.-x„^„_i): Xij e ]R2 foj. i^jj^ (7) 

where we call all the variables x^ j for a particular index 
i the replicas of the original variable x^ . Every configura- 
tion (xi, . . . x„) G can be identified with a configura- 
tion (xi^2, ■ • -Xn^n-i) G 17, whcrcin x^.j = x.; for all 
through a simple linear map A. Enough redundant vari- 
ables have been introduced to the configuration space so 
that each pairwise exclusion constraint can now be writ- 
ten in terms of a private set of variables, disjoint from 
the private variables of other constraints: 



D''^ = {(xi^2,---x„,„_i) e 17: ||x, 



^,3 



>2}. (8) 



* = {(xi,X2):xi,X2eR2}. 



(1) 



The intersection, DC J7, of all of the pairwise exclusion 
constraints, which we will call the "divide" constraint, is 
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(a) 








FIG. 1: An illustration of the D — C scheme in the case 
of packing three disks into a square box. In the case of two 
overlapping disks (a), the exclusion constraint is simple in the 
sense that a projection to the constraint can be performed ef- 
ficiently. The projection is given by ((3| and yields the config- 
uration (b). In the case of three disks (c), there is no similarly 
efficient projection method to the constraint that no overlaps 
occur. In the D — C scheme, each disk is represented by two 
replicas (d), which together make three independent replica 
pairs (e). The exclusion constraint, now also called the "di- 
vide" constraint, is modified so that only replica pairs are 
prohibited from overlapping, and any other overlaps are al- 
lowed. Thus, the projection to the exclusion constraint can be 
performed independently on each replica pair as in the case of 
two disks (f). A second constraint, the "concur" constraint, 
requires all replicas representing a single disk to coincide and 
requires the disk to lie within the confinement box. The re- 
sult of projecting the configuration (d) to this constraint is the 
configuration (g). In order to search for a configuration sat- 
isfying both constraints, we do not alternately project from 
one constraint to the other, but instead use the difference 
map (|12p to evolve the configuration. The result of a success- 
ful search is a configuration (h) satisfying both constraints, 
which by construction corresponds to a solution of the prob- 
lem. 



now also simple, since the projection can be performed 
independently on each set of private variables (Figure [ij. 

The map A : ^ Q, from the original configuration 
space (the physical configuration space) to the new one 
(the formal configuration space) is not surjective, and so 
a general point in the formal configuration space does not 
correspond to a vahd physical configuration. The "con- 
cur" constraint C = A(^') is given by the range of the 
subset of ri that does correspond to vaUd configurations. 
That is, the constraint requires redundant specifications 
of an original variable to concur in regard to its value. 
Since A is linear, C is also a simple constraint. 

Another constraint that must usually be addressed in 
packing problems with a finite number of particles is the 
confinement constraint. In most cases the particles, or 
their centers, are confined to lie in some subset M of 
space, where M can be either some region of finite vol- 
ume, or a compact manifold (as in the case of spherical 
codes). As a subset of the original configuration space, 
the confinement constraint is written as 

i^conf = {(xi, . . .x„) e^:^,eM for ah i} C (9) 

We can incorporate this constraint into the "concur" con- 
straint, C, by modifying it to be the image .A(_/^conf) 
stead of the entire range of A. In our example, this would 
give the constraint 

C {(xi.2, . . .x„.„_i) e fJ: = Xi e M for aU i, j}. 

(10) 

Since A is linear, the projection to C = A{Kconi) can 
be decomposed into a projection to A{'^) followed by 
a projection to A{Kconf)- The first step is performed by 
taking the average position of all the replicas of each disk. 
The second step is performed by projecting this average 
position to M (see Figured]). In general, this two-step 
projection method is valid for handling the constraints 
in the physical configuration space that are simple at 
the outset and do not require the introduction of new 
variables. 

The result of the above construction is that a configu- 
ration in n satisfies the "divide" and "concur" constraints 
simultaneously if and only if it corresponds to a config- 
uration in ^E* which satisfies all the exclusion constraints 
and the confinement constraint; that is, it corresponds to 
a solution of the packing problem under consideration. 

In the following sections we modify the above simple 
construction so as to generalize the method in two major 
ways. The first generalization is to packings of infinite 
regions, instead of only finite ones. Specifically, we allow 
for periodic packings with an arbitrary unit cell. This 
is achieved by generalizing the idea of replicas of parti- 
cles to include also their periodic images. When the unit 
cell vectors are included in the original set of parameters, 
the map from the original parameter space to the space of 
replica configurations is still linear, though a little more 
elaborate. Additionally, the confinement constraint of 
finite packings is replaced in the case of periodic pack- 
ings by a constraint on the unit cell volume, ensuring a 
specified density. 



5 



The second generalization is to packings of non- 
spherical particles, specifically convex polytopes. This 
is achieved by representing each particle not only by the 
position of its centroid, but by the positions of all its ver- 
tices. A new constraint, the rigidity constraint, is added 
to ensure that the particle is not deformed in the solution. 
Despite the mathematical complications that arise from 
these two generalizations, the conceptual framework is 
identical to the above example, and the constructions in 
the following sections will draw attention to the analogy 
with the construction presented above. 



B. The difference map 

Given a problem formulated as the task of finding a 
configuration x G C H D, simultaneously satisfying the 
constraints C,D C fj, we wish to use the availability of 
efficient methods for computing the projections ttc and 
TT^i to the constraints in order to set up an iterated map 
to search through the configuration space for a solution. 
Naive schemes, such as the alternating projections map 
X H- 5- 7r£i(7rc(x)), suffer from the problem of stagnation at 
near solutions (local minima of the distance between the 
two constraints). The difference map, a slightly more so- 
phisticated scheme, is designed to provide efficient search 
dynamics while avoiding the traps of local minima (isj . 

The difference map (DM) can be written in terms of 
the projections nc and ttd and one parameter /3: 

DM:fl^n (11) 
XK^x + /?[7rz5(/c(x))-^c(/D(x))], (12) 

where 

/d(x)= (l-^) 7rz3(x)-hix, 



/c(x)= + 7rc(x)--x. 

In this paper we use only /3 — 1. A difference map 
search proceeds by starting from a random initial config- 
uration xq and iteratively applying the difference map: 
Xi = DM(xi_i) fl^. When the map reaches a fixed point 
x/p, a solution is obtained by 

Xsoi = TTC (/d(x/p)) = TTD (/c(x/p)) . (13) 

Notice that the ability to obtain a solution from any fixed 
point of the map, due to the cancelation of the two brack- 
eted terms in , relies on the definition of the problem 
in terms of only two simple constraints. For a given it- 
erate Xi, the terms ttq (foi^i)) and ttd (fci^i)) provide 
two estimates of the solution, each satisfying one of the 
two constraints. We call these respectively the C- and 



D-estimates of the solution at the ith iteration. The dis- 
tance between the two estimates is the error e and the 
search terminates when the error converges to zero. 

To summarize, a simple difference map solver for con- 
tinuous constraints would consist of the following simple 
steps: 

1. Initialize the iterate x to a random configuration. 

2. Compute the two estimates of the solution xc 
T^c (/r>(x)) and x/, ^ ttd (/c(x)). 

3. Compute the error e ^ ||xc — X£)||. If it is be- 
low a predefined convergence threshold, the search 
terminates, and the solution is given by xp w x^. 

4. Advance the iterate x -s— x -I- /3(xd — xc). Start 
the next iteration at Step 2. 

III. CONSTRAINTS 

A. Periodic sphere packing and kissing 

A periodic packing of equal-sized spheres (radius r) in 
d dimensions can be generated by the action of a lattice 
A on a set of p primitive spheres. Let P be the set of 
centers of the primitive spheres. We define a generating 
matrix of the packing as a (d + p) x d matrix M whose 
first d rows are a set of generators of A and whose re- 
maining p rows are the vectors in the set P. Combining 
these quite different sets of configuration variables into 
a single matrix serves to remind us that at the highest 
level of our search algorithm both sets are treated in a 
uniform manner by the projection operators. The de- 
tailed constraints, of course, distinguish among the two 
parts of M, which we denote Mq (lattice generators) and 
Ml (primitive sphere centers). The set of all the centers 
of spheres in the packing is then the Minkowski sum 

A + P = {boMo+y: bo eZ^ yeP} (14) 
= {bM: b e Z'^iSEp}, 

where Ep is the set of coordinate-permutations of the p- 
dimensional vector (1,0, 0, . . . , 0). The space M(^+p)><'^ of 
generating matrices takes the role of the physical config- 
uration space 

A matrix M generates a valid packing if the centers 
of any two sphere of the packing, biM and b2M, are 
separated at least by a distance of 2r when bi 7^ b2. 
Each choice of bi and b2 generates a constraint on the 
matrix M 

llbiM-bsMII > 2r, (15) 

which we call an exclusion constraint. Note that there are 
infinitely many independent exclusion constraints (con- 
straints with bi — b2 = h'^ — h'2 are not independent). 
However, for any non-degenerate matrix M only finitely 
many independent exclusion constraints are violated or 
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are even remotely close to being violated. In practice, 
only those constraints need be tested in our computa- 
tions. We call those constraints the relevant exclusion 
constraints (let there be n of them), and we define a 
2n X {d + p) matrix A whose rows a2i-i and a.2i are 
the vectors bi and b2 related to the ith relevant ex- 
clusion constraint. We discuss below how the relevant 
constraints are identified. 

The linear map ^ : M X = AM is a map from the 
physical configuration space 'J to a larger-dimensional 
space, — M^"^*^, which we use as the formal configura- 
tion space. As before, since the map A is not surjective, 
only a subset (a linear subspace, in fact) of formal config- 
urations have a corresponding generating matrix in the 
physical configuration space. The choice of constraints 
below will guarantee that solutions belong to this subset. 
The size of the configuration space grows as the number 
of relevant independent exclusion constraints, which is 
the number of independent near neighbor pairs for which 
overlap needs to be actively avoided. Notice that each 
relevant exclusion constraint can now be written in terms 
of a private set of variables. Specifically, each row of the 
matrix X corresponds to the position of one particle, and 
the ith relevant exclusion constraint is given by 



A = {Xer!: I |x2.-i -X2.il > 2r}. 



(16) 



The intersection of all the relevant exclusion constraints 
forms our "divide" constraint, 

D = {X e f2: ||x2i-i -X2,|| > 2r for i = 1,. . .n}. (17) 

Each set of private variables associated with one exclu- 
sion constraint is composed of the coordinates of replicas 
of two particles, and we call these two replicas a replica 
pair. 



As mentioned in Section III Al the confinement con- 
straint of finite packing problems is replaced in the case 
of periodic packings with a constraint on the density of 
the packing. The density of a packing generated by a 
matrix M is given by the density of the unit cell, whose 
volume is | det Mq | and which contains p particles of vol- 
ume Vi. 



I det Mn 



(18) 



Therefore, if we wish to find a packing of density </> > 
'/'target, the density constraint on the generating matrix 
will be 



density 



{Me IdetMol < V^argct}, (19) 



where Vtargct = pVi/(j>t argct ■ 

As in the example of Section III A| since the map A 
is not surjective, a general element X e of the formal 
configuration space does not correspond to a well-defined 
physical configuration. We therefore impose a constraint 
that requires X to lie in the range of A. In the context of 
the PDC construction we call this the lattice constraint 



because it requires different periodic images of a primitive 
particle to lie on the points of a lattice, and requires that 
lattice to be the same for all primitive particles (up to 
translation). Again, as in Section FlI Al we combine the 
lattice constraint with the density constraint to form the 
"concur" constraint: 



C* — ^(-ft^dcnsity) 

={X ^ AMefl: 



(20) 



IdetMnI < Vt, 



With these definitions of the constraint sets, X = 
AM G CnD if and only if M generates a periodic packing 
of density 4> > ^target- The action of the projections ttd 
and TTc to the two constraints is illustrated in Figure [2] 
and Sections IIV Al and IIV BI discuss how the projections 
are computed efficiently. 

The basic operations of the search — projections — 
depend directly on the metric defined on the formal con- 
figuration space. Therefore, the choice of metric affects 
both the complexity of implementing the projection and 
the search dynamics. The simplest choice for the metric 
is the distance induced from the Frobenius (Euclidean) 



|Xi - X2III, = trace ((Xi - X2)(Xi - X2)' 



(21) 



This choice of metric amounts to giving all replicas of a 
particle equal weight in influencing its consensus position 
in the "concur" projection. We can use a slightly different 
Euclidean metric, given by 



|Xi 



X2IIW 



trace (W(Xi - X2)(Xi - X^ 



(22) 



where W is a diagonal matrix whose diagonal elements 
Wi are the metric weights of different replicas. Per- 
formance is greatly enhanced by adjusting the metric 
weights throughout the search to afford greater weight 
to replica pairs that continually violate their constraints 
and smaller weight to replica pairs that are in low risk 
of violating their constraints '19]. Note that removing a 
constraint from the list of relevant constraints (i.e., re- 
moving the corresponding pair of rows from A and X) is 
equivalent to setting the metric weight of its replicas to 
zero. Therefore, in the course of the search we not only 
adjust the weights Wi of replica pairs, but also add and 
remove replica pairs. The details of how these changes 
are applied systematically are given in Section [IV CI 

This constraint formulation of the periodic sphere 
packing problem (finding a periodic packing with density 
^target) can be straightforwardly modified to describe in- 
stead the periodic kissing number problem (finding a pe- 
riodic packing with average coordination number Ttaigct)- 
First, the "divide" constraint is modified so that each 
replica pair must still be separated by a distance of at 
least 2r, but at least prtarget replica pairs must be sep- 
arated by a distance of exactly 2r. Second, the condi- 
tion on the volume of the unit cell is dropped from the 
"concur" constraint. Projections to these modified con- 
straints are also given in Section IIVI 
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FIG. 2: An illustration of the "divide" and "concur" pro- 
jections in the two-dimensional sphere packing problem with 
p = 3. (a) A hypothetical configuration of six replica pairs in- 
volving the primitive disk A. Disks with the same letter mark- 
ing their centers are replicas of the same primitive disk (as for 
disk A) or of its lattice translates (as for disks B and C). One 
replica pair, violating its exclusion constraint, is emphasized, 
(b) The output of the "concur" projection: the closest config- 
uration to (a) such that all replicas of a particular primitive 
disk lie on top of each other, or a lattice translation apart (ar- 
rows), and such that those lattice translations define a lattice 
with a sufficiently small unit cell volume. This projection is 
a modification of the "concur" projection depicted in Figure 
[TJl,g. (c) The output of the "divide" projection: the clos- 
est configuration to (a) such that no replica pair violates its 
exclusion constraint. This is identical to the "divide" projec- 
tion depicted in[TJi-f. Detail: (d) the emphasized replica pair 
before the "divide" projection (thin-outline disks) and after 
(thick-outline disks) isolated for clarity. 



B. Convex polytope packing 

The symmetry of the spherical particle allows its con- 
figuration to be described solely by the position of its cen- 
ter. In the case of a general convex particle, the variables 
of the configuration space need to include information 
also about the orientation of the particle. One possible 
description of the particle assigns variables separately to 
the position of its centroid and to the description of the 
rotation about the centroid (e.g., a rotation matrix or a 
quaternion). In this paper, however, we find it more con- 
venient to describe convex polytopes by reference to the 
positions of their vertices. Therefore, a polytope with 



V vertices is represented hy a v x d vertex matrix and 
is given by the convex hull of these vertices. Although 
the configuration of a single particle is no longer rep- 
resented by a single vector but by a matrix composed 
of V vectors, it is convenient to treat these matrices as 
vectors, which we typeset as bold-face upper-case Latin 
letters (e.g., X for the vertex matrix of the polytope 
K = convX = convjxi : i = 1, . . . u}), and to construct 
matrices whose rows are such vectors. A translation by 
t of a polytope convX is given by conv(X -I- c-^t), where 
c'^ is a column vector of unit elements and c-^t is the 
translation matrix corresponding to the translation vec- 
tor t. Similarly, a rotation is given by conv(XR), where 
R is a d X d orthogonal matrix. 

A periodic packing is again generated by the action of 
a lattice A on a set of p primitive polytopes whose vertex 
matrices form the set P. The set of all vertex matrices 
of polytopes in the packing is the Minkowski sum 

c^A + P = {boMo + Y : bo e Z'', Y e P} (23) 
= {bM: b e Z'^eEp}, 

where M is a generating matrix of the packing, whose first 
d rows (comprising Mq) are translation matrices generat- 
ing A, and whose remaining p rows (comprising Mi) are 
the vertex matrices of the set P. The space of generating 
matrices ^' = ]R.('^+p)^('"^'') is the physical configuration 
space. 

Each exclusion constraint between two particles of 
the packing requires the convex hulls conv(biM) and 
conv(b2M) not to overlap for any bi ^ h-z- To con- 
struct the formal configuration space we again form one 
replica pair for the particles involved in each relevant ex- 
clusion constraint, which gives fl = I|j2nx(uxd)^ rjij^^ map 
A from physical configurations to formal configurations 
is given by the matrix A whose rows a2i_i and a.2i are 
the vectors bi and b2 related to the ith relevant exclu- 
sion constraint. The "divide" constraint is given by the 
intersection of all the relevant exclusion constraints, each 
expressed in terms of its private replica pair: 

D = {XeQ: int(convX2i-i n convX2i) = (24) 

for z = 1, 2, . . . n}. 

In addition to the lattice constraint and the density 
constraints, which combine in Section IIII Al to form the 
"concur" constraint, in the case at hand we must include 
a third constraint, the rigidity constraint. The primitive 
particles of a packing generated by a general matrix M 
are only constrained in their number of vertices, not in 
the arrangement of those vertices. However, we are inter- 
ested only in packing where all the particles are congruent 
with a given shape, and so we impose the constraint on 
M that the vertices of its primitive particles are obtained 
from the vertices of the given particle by a rigid motion: 

i^rigidity = {M e * : Y = Y(o)R, + c^t, (25) 
for all p rows Y of Mi}, 
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"divide" constraint 


"concur" constraint 


sphere 
packing 


X2i-i — X2i > 2r for 
all n replica pairs 


X = AM 

M G A'density 


kissing 
number 


||x2j_i - X2i|| > 2r for 
all n replica pairs and 
= 2r for prtargct pairs 


X = AM 


polytope 
packing 


convex hulls of X2J-1 
and X2i+i non-overlap- 
ping for all n replica 
pairs 


X = AM 

M £ /Cdcnaity 

and M G ifrigidity 



TABLE I: A summary of the D — C constraints for peri- 
odic sphere packing, the average kissing number problem, and 
polytope packing. The "divide" constraint encompasses the 
relevant exclusion constraints, while the "concur" constraint 
encompasses, where applicable, the density, rigidity, and lat- 
tice constraints. 

where Y^"^ is the vertex matrix of the given particle. The 
"concur" constraint C = ^(ii'dcnsity H i^rigidity) is given 
by combining the density and rigidity constraints on the 
generating matrix with the lattice constraint. The result 
of constructing the "divide" and "concur" constraints is 
that a formal configuration satisfies both of them if and 
only if it corresponds to a generating matrix in the physi- 
cal configuration space which yields a packing of the given 
particle with the desired density. Table |T] summarizes the 
D — C constraints for the three problems discussed and 
Section |IV] describes in detail the projections to these 
constraints. 



IV. IMPLEMENTATION 
A. "Divide" projections 

1. Sphere packing and kissing 

In order to implement an iterated difference map 
search, whose iterations are given by (|12p . we must im- 
plement efHcient projections to the "divide" and "con- 
cur" constraints. These implementations are the subject 
of Sections IIV Al and IIVBI In the course of the search, 
considerations of efficiency require certain changes to the 
formal configuration space - specifically adding and re- 
moving replica pairs, changing metric weights, and lat- 
tice reduction. In section flV CI we discuss when and how 
these changes are applied. 

In the case of sphere packing, the "divide" constraint 
simply requires that the centers of the two spheres com- 
prising each replica pair be a certain distance apart. This 
is obtained by applying equation ([3]) to each replica pair. 
Note that the "divide" projection acts independently on 
each replica pair, and since the metric weight of all vari- 
ables specifying one replica pair are equal, the metric 
weights have no influence on this projection. The action 
of this projection is illustrated in Figure [5] For the kiss- 
ing number problem, the first case of ([3]) is also used if 




FIG. 3: An illustration of the polytope exclusion and rigidity 
constraint projections for the case of regular pentagons, (a) 
The pentagons are non-overlapping, as demonstrated by the 
existence of a separating axis that passes through one vertex 
of each pentagon (red dots), (b) A hypothetical situation 
of overlapping pentagons. No axis, and particularly no axis 
passing through one vertex of each pentagon, separates the 
two sets of vertices, (c) For the input pentagons in (b), the 
subset S — T (red dots) that minimizes S{S)^, the sum of 
squared-distances to the least-squares axis (red line) while 
satisfying that the latter separates the remaining vertices, (d) 
Another choice of S that yields a valid separating axis, but a 
larger sum of squared-distances. (e) A choice of S that yields 
an axis that fails to separate the remaining vertices, (f ) Using 
T found in (c), the output of the projection to the exclusion 
constraint is determined by moving the points of T onto their 
least-squares axis, (g) The output (solid line) of the rigidity 
projection for the input pentagons (dashed) in (b). 



the replica pair is one of the prtarget closest replica pairs. 



2. Convex polytope packing 

Although identifying and resolving overlaps between 
two spheres is straightforward, the same task is more 
challenging in the case of other convex objects, where 
more degrees of freedom come into play. The literature 
on the topic of detecting overlaps (collisions) between 
polyhedral solids is extensive (driven in part by appli- 
cations in computer graphics), and many efficient tech- 
niques exist for checking whether two convex polyhedra, 
convXi and convX2, overlap (see e.g., [H, Hj])- In our 
case, as we are interested in computing the projection 
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to the exclusion constraint, we also need to determine 
the distance-minimizing resolution of the overlap. That 
is, we must find the smallest displacement of the vertices 
such that the new polyhedra do not overlap. As far as we 
have been able to determine, there is not an established, 
efficient computational method developed for this specific 
problem. The method we provide here is efficient enough 
for the purpose of packing polyhedra with a small num- 
ber of vertices, but the computation time required grows 
exponentially with the number of vertices. A more effi- 
cient resolution method for particles with more vertices 
and for smooth particles is currently in development. 

The method relies on the separating plane theorem: 
the convex hulls of two sets of vertices in K'' do not over- 
lap if and only if there is a (d — l)-dimensional plane 
that separates the two sets, so that each is contained in 
a different half-space. The theorem can be made even 
stronger by specifying that the separating plane can al- 
ways be chosen to contain d vertices from the given sets, 
including at least one from each. Therefore, one can 
check whether two polytopes overlap by checking whether 
they are separated by any of the planes defined by any 
such subset of vertices (Figure |3^-b) . If the polytopes 
are non-overlapping, the resolution leaves the vertices un- 
changed. 

If the polytopes overlap, we must find the smallest dis- 
placement of their vertices that resolves the overlap. In 
the resolved configuration there is a separating plane, 
Vsp = {r G R'^ : fisp -r = /isp}, that separates the two sets 
of vertices. As a consequence of distance minimization, 
the only vertices moved in the course of the resolution are 
the vertices which lie on Vsp in the resolved configuration. 
Let T and T' be the pre- and post-resolution positions, 
respectively, of those vertices that are displaced during 
the resolution. Therefore, T' is the set of points in Vsp 
closest to the points of T: 



r = {r+ {hsp - r • nsp)nsp : r e T} C V 



' sp ■ 



(26) 



The squared norm of the resolution displacement is 



(27) 



For any set of points 5* there is at least one plane V — 
{r e M"* : n • r = ft.} that minimizes the sum of squared 
distances 



res 



(n • r - hf 



(28) 



We call such a plane a least-squares plane of S. The 
separating plane of the resolved configuration is always 
a least-squares plane of T. If this were not the case, a 
small tilting of the separating plane towards such a least- 
squares plane (with a corresponding movement of the 
points in T') would result in a resolution by a smaller dis- 
placement. In order to resolve an overlap between poly- 
topes conv Xi and conv X2 , we therefore have to solve a 
discrete problem: among all subsets S' of ATi U X2 with 



a least-squares plane separating the remaining vertices 
Xi \ S from X2 \ S, find the one with the minimal sum of 
squared distances (pSj) . This is the set T (Figure [5t-f). 

The least-squares plane Vs of a set S is determined 
by minimizing the sum of squared distances (j28l) . Note 
that for a fixed normal direction n, the value of h that 
minimizes the sum is h — fi -r, where f = X^res ^/\'^\ 
the centroid of S. Therefore, we wish to minimize 



5^[n.(r-r)]2 = n 



res 



^(r-r)'^(r-r) 



.res 



n , 



(29) 



the minimum of which is equal to the smallest eigenvalue 
of the symmetric matrix J2res(^ "^^^ mini- 

mum is realized when n is the corresponding eigenvector. 
Degenerate cases with equal lowest eigenvalues occur, but 
they do not pose a problem: whenever an optimal sepa- 
rating plane occurs as a degenerate least-squares plane of 
some set S, its degeneracy implies that there is a least- 
squares plane of S which also includes an extra vertex; 
this plane will be equally optimal and will occur as a 
less degenerate least-squares plane of a superset S' 3 S. 
Therefore, the optimal least-squares plane always occurs 
as a non-degenerate least-squares plane of a set S. 

To summarize, the overlap detection and resolution al- 
gorithm consists of three steps (illustrated in Figure [5^- 
f): 

1. Consider all subsets S C Xi U X2 of size |5| — d 
with at least one point from each polytope. Let 
V — {r Kf^ : h ■ r — h} he a, plane that includes 
S. For each S let 



AliS)^ ^(n.x-ft)2+ ^(n.x-ft)2, (30) 

nx>/i nx<^ 

AUS)^ ^(n.x-ft)2+ ^(n.x-ft)2, (31) 

(32) 



nx</i 



xeX2 

n-x>/i 



A2(5)=min(A^(5),A2_(5)), 
and let 



=minA^(S'). 

s 



(33) 



A^ provides a measure for the interpenetration of 
the two polytopes. If A^ — 0, then a separating 
plane exists, the input polytopes do not overlap, 
and the algorithm ends here by returning the orig- 
inal vertex positions Xi and X2- If A^ > 0, the 
polytopes overlap and the algorithm continues to 
Step 2. 

2. Consider all subsets 5 C Xi U of size \S\ > d 
with at least one point from each polytope. Let 
= {r g R'' : n • r = /i} be a least-squares plane 
of S*. If the plane separates the vertex sets with 
the points of S removed — Xi\ S and X2\ S — 
let S'^{S) be the sum of squared-distances from S 



10 



to the plane. Otherwise, let S^{S) = oo. Among 
the subsets S considered, let T be the subset that 



minimizes S (S) and Vr — {r 'EM. : n^-r 



hx} be 



its associated least-squares plane. As (5^(5) < oo 
if S contains all vertices, the minimum is always 
finite. Continue to Step 3. 

3. The sets of vertices returned are given by X[ and 
X2 , wherein x' G X[ U X2 is given by 




if x^r 

X • XiT^iiT if X G T, 



(34) 



where x G Xi U X2 is the corresponding original 
vertex position. 

The projection 7r£)(X) to the "divide" constraint (l24l) . 
of an input matrix X comprised of pairs of vertex matrices 
X2i_i and X2i, is then achieved by applying the above 
algorithm independently to all i = 1, ... n pairs. 



B. "Concur" projections 

1. Lattice constraint 

All the "concur" constraint sets described in this paper 
are of the form 



C = A{K) = {X = kMeVt-.Ue K) 



(35) 



where A is constant, and M is variable, but must satisfy 
a constraint M € i^. The projection then is given by 

TTc : X X' = AM, (36) 

where M realizes the minimum over K of the distance 

\\X-X'\\^ = trace (W(X- AM) (X-AM)^) . (37) 

Absent any constraints on M (as for example in the 
"concur" constraint for the kissing number problem, 
where if = Vl/), the solution would be given by 



M = (A^WA)~iA^WX. 



(38) 



This can easily be seen by writing M = M + (5M, which 
gives 



\X- X 



= trace (W(X - AM)(X - AM)^) 
= c + trace(WA 5M 5M^ A^) 
= c + trace(W'(5M 5M^), (39) 



where W' = A-^WA and the constant term c does not de- 
pend on (5M. The second term is non-negative, and when 
M is unconstrained, (15^1) is minimized by letting M = M. 
Additionally, we have just reduced the constrained case 
to the problem of finding M G if that minimizes the cost 
function 



/(M) = trace (W'(M - M)(M - M)^) 



(40) 



This projection strategy parallels the two-step strat- 
egy used in Section fll Al First, the formal configuration 
X is projected to the range A{'^) of the physical config- 
uration space, giving AM. Then, the projection of M to 
the additional constraint K is performed in the physical 
configuration space using the metric induced on its im- 
age in the formal configuration space. Below, we solve 
the second step of this projection problem for various 
constraints K. 



2. Density constraint 

In the "concur" constraint for the sphere packing prob- 
lem, the only constraint on the generating matrix is the 
density constraint. The set of generating matrices M sat- 
isfying the density constraint is 



density 



= {Af: IdetMol < Vtargct}, 



(41) 



where Mq is the generating matrix of the lattice and is 
given by the first d rows of M. If | det Mo| < Vtargct, then 
the projection to the constraint (the choice of M that 
minimizes the cost function (|40|)) is trivially M = M. 
Otherwise, since Mi is unconstrained, we can minimize 
with respect to Mi for a given Mq. This yields Mi = 



Ml 



W'li Wio(Mo 



Mo), where Wjj are the block- 
elements of W acting on M/ to the left and on Mj to 
the right. Thus, the cost function for Mq is simply 



/(Afo) = trace (W"(Mo - Mo)(Mo - Mo)' 



(42) 



where W" = W^q - Wf^iW^i^W; 
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The projection becomes easier to analyze in terms of 
the matrix L — (W")^/^Mo. The cost function then takes 
the form of the simple Frobenius distance 

/(L) = trace ((L-L)(L-L)^) , (43) 

and the density constraint is still in the form 



|detL| < V4g3t, 
Idetl^'T/". 



(44) 



where V4i.got = Kaiget/| det Since the absolute 

value of the determinant of L is given by the product 
of its singular values, the solution to this minimization 
problem is given by a matrix L = USV with the same 
(right and left) singular vectors as the matrix L = UEV, 
but different singular values. The cost function expressed 
in terms of the singular values Ui and a7 of, respectively, 
L and L takes the form 



/(E)=^(a.-a 



(45) 



We numerically minimize this quadratic function subject 
to the density constraint (|44l) . Through back substitu- 
tion we then have the matrix M that minimizes (1371) and 
TTc{X) = X' = AM. 



11 



Rigidity constraint 



C. Formal configuration space maintenance 



In the "concur" constraint for the polytope packing 
problem, an additional constraint on the generating ma- 
trix M is that the primitive polytopes that make up Mi 
are congruent with a given polytope. The generating 
matrix is then constrained to the set 



K = K. 



density 



rigidity 



(46) 



where 



i^dcnsity = {M: I detMol < ^target}, 

i^rigidity = {M: Y = Y(0)R,+c^ti 
for all p rows Y of Mi}. 



To calculate the projection ttcOC), the cost function (|40)) 
must be minimized over K. However, since the off- 
diagonal block Wqi couples the lattice parameters Mg 
to the primitive particle parameters Mi, this minimiza- 
tion is complicated. Instead of exact minimization, we 
employ a two-step heuristic method, which results in an 
approximate projection. 

In the first step, we calculate the matrix M' g -K^donsity 
that minimizes the cost function, as in Section IIVB 21 
Then, in the second step, we calculate the matrix M G K 
by applying to each row Y of M'l the smallest change so 
that it becomes a vertex matrix of a polytope congruent 
with the reference polytope. The second step is achieved 
by finding the rigid motion applied to the reference poly- 
tope which brings its vertices as close as possible to the 
vertices of Y as measured by the sum of squared distances 
(Figure[3^). The problem of finding the rigid motion that 
brings one given list of points closest to another given list, 
sometimes known as the problem of absolute orientation, 
occurs frequently in a variety of fields (e.g., in calculat- 
ing RMSD between two conformations of a biomolecule) 
and several efficient methods for its solution have been 
developed (see [Ull^). 

The output of the approximate projection is then given 
by X' = nc{X) = AM « 7rc(X). As X' e C, the approx- 
imate projection gives a configuration in the constraint 
set, but might not give the closest one to the input config- 
uration. We justify the use of the approximate projection 
by noting that it is an exact projection if the off-diagonal 
block Wqi is zero. A non-zero off-diagonal block is the 
result of correlations in the relevant exclusion constraint 
vectors b between the coefficients of lattice translations 
and the coefficients of primitive particle vertex positions. 
We expect these coefficients to give uncorrelated contri- 
butions and to add up to small off-diagonal elements due 
to random cancellations. Indeed, we find that the off- 
diagonal block is small in comparison with the diagonal 
blocks, and we expect our heuristic to yield a good ap- 
proximate projection. 



In our discussion of the choice of metric in Section Hill 
we discussed the ideas of dynamically readjusting the 
metric (through the weights Wi of the various replicas) 
and of removing and adding replicas (removing replicas is 
formally equivalent to setting their weight to zero). The 
latter is necessary for implementation reasons: there are 
infinitely many independent exclusion constraints (and 
therefore replicas), but we can only represent a finite 
number of replicas in our implementation. As the set of 
relevant constraints changes over the course of the search, 
we must remove and add replicas. Our criterion for which 
replicas to represent is based on the difference map's cur- 
rent "concur" estimate: we include a replica pair for each 
pair of particles whose centroids in the "concur" estimate 
are closer than some cut-off distance. Using the gener- 
ating matrix obtained in the "concur" projection we can 
easily find all such pairs using the method of Agrell et al. 
p?! ]. The cut-off distance is chosen so that at least all 
replicas that might be in risk of overlap are represented. 

The problem of implementation is not the only reason 
we wish to limit the number of replicas we represent. A 
proliferation of unnecessary replicas has the adverse effect 
of attenuating the information obtained from the "con- 
cur" projection by diluting the influence of more critical 
replicas. We observe that such replica proliferation could 
result not only in a slower search, but also in an increased 
tendency to become trapped in local optima. Limiting 
the number of replicas is one way to avoid this effect, 
but we find it useful to further amplify the information 
from critical constraints by giving them greater weights 
[l9| . We perform the weight adjustments adiabatically, 
that is, slowly over the course of many iterations, by up- 
dating the weights of each replica pair according to the 
rule 



TWi 



Wi 



1 



(47) 



where w-(Xc) is a function that assigns replicas weights 
based on their configuration in the "concur" estimate, 
and T is a relaxation time for the replica weights in units 
of iterations. 

In the sphere packing problem (in d dimensions, with 
unit spheres), we choose the weight function to be 




l|x,||^) 

2 _ g-j-2-d/2 



if 
if 



< 2 
> 2, 



(48) 



with a 20. The dimensional dependence is chosen so 
that under the assumption of uniform density, the to- 
tal weight from replicas over a certain distance follows a 
dimension-independent power law. In the polytope pack- 
ing problem, we similarly use 




if the polytopes overlap 
-4r|„)-2 if not, 

(49) 
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with a « 10, where rm is the inradius of the polytope, 
ri is the centroid-centroid distance of the polytopes, and 
Af is the measure of the overlap between the polytopes 
defined in ([55)1 . 

In addition to the maintenance of rephcas, which is 
performed after every iteration of the difference map, 
we also periodicahy perform a lattice reduction using 
the LLL algorithm ^2^]. The lattice generated by Mq is 
re-represented using the LLL-reduced generating matrix 
Mq = GqMo, where Go is a unimodular integer matrix. 
Additionally, all primitive particles whose centroids are 
outside of the unit ceU given by A^a^ : — 1/2 < < 
1/2} are re-represented by their lattice-translate in that 
cell. In summary, the new packing generating matrix M' 
is given by 

M' = GM = J ^ M, (50) 

where Gi gives the lattice translations to be applied to 
the primitive particles. Since the actual positions of the 
particles, as represented in the matrix X = AM, should 
be unchanged, the lattice reduction must also be applied 
to the nominally constant matrix A (A — > A' = AG~^). 



V. RESULTS 
A. Sphere packing 

Using the PDC scheme described in the previous sec- 
tions we perform a de novo search for the densest lattice 
{p = 1) sphere packings in dimensions 2 — 14. The PDC 
search, starting from random initial configurations, was 
able to reproduce the densest packing lattices known for 
all cases, and the results of the search are summarized in 
Table [Til For dimensions 2 — 8 the lattices are known to 
be optimal, and for dimensions 9 — 14 these results are, 
to our knowledge, the first numerical evidence from a de 
novo search that the known lattices are optimal. 

Note that the number of replicas is determined by the 
number of near neighbors of each sphere, which rises 
rapidly with the number of dimensions. This rise causes 
an increased computational storage cost per physical de- 
gree of freedom in a PDC search, compared to a con- 
stant storage cost per physical degree of freedom in a 
method involving a local search in the physical configura- 
tion space. However, this rise need not affect the scaling 
of CPU costs, since both search methods need necessarily 
check a comparable number of particle pairs for possible 
overlaps. 

In dimensions d = 10,11,13 there are known non- 
lattice packings with p — 40, 72, 144 respectively that 
are denser than the densest known lattices [11] . In up to 
11 dimensions, we searched for non-lattice packings with 
as many as p — 12 primitive spheres, but the searches 
did not produce packings denser than the lattice pack- 
ings. For a density target matching the lattice density. 



d -^densest 


6 

T^dcnscst 


{Niter) 


(n) 


titer success rate 


2 A2 


0.90690 


42 


11 


0.1ms 100/100 


3 D3 


0.74047 


230 


38 


0.2ms 100/100 


4 D4 


0.61685 


191 


127 


0.4ms 100/100 


5 D5 


0.46526 


308 


323 


1ms 100/100 


6 Ee 


0.37295 


173 


977 


2ms 100/100 


7 Et 


0.29530 


217 


2740 


5ms 96/100 


8 Es 


0.25367 


99 


8528 


20ms 96/100 


9 Ag 


0.14577 


161 


16314 


30ms 85/100 


10 Alo 


0.092021 


394 


31433 


70ms 47/100 


11 Kii 


0.060432 


421 


68722 


0.3s 54/100 


12 K12 


0.049454 


397 


204321 


0.9s 55/100 


13 Ki3 


0.029208 


577 


430796 


2s 25/100 


14 Ai4 


0.021624 


1652 


1007250 


6s 4/10 



TABLE II: Results of PDC searches for dense lattice packing 
in dimensions d = 2, ... 14. For each dimension, 100 runs from 
random initial conditions were performed with the density 
target ^target densest' dcusity of the densest known 
lattice Adensest The runs were limited to 5000 iterations, 
and the number of converged runs is quoted in the right- 
most column. For dimensions 10 and above, each run was 
first allowed to converge at a density target of O.80densest and 
then continued with the final target. The mean number of 
difference map iterations in converged runs was (A^'iter), and 
the mean number of relevant exclusion constraint used was 
(n). Each iteration took an average runtime of titer on a 
single 3 GHz CPU. In d = 14 only 10 runs were performed 
with three intermediate targets. 



the searches reproduced the lattice packing, suggesting 
that the lattice packing in these dimensions is the opti- 
mal packing with a small number of spheres in the unit 
cell. 



B. Kissing number 

For the kissing number problem, PDC searches were 
able to reproduce the best known lattice kissing arrange- 
ments in dimensions 2 — 11. In dimensions 2 — 9, the re- 
sult is known to be optimal, and for dimensions 10 and 
11, we are not aware of previous numerical evidence for 
their optimality. Table IIIII summarizes the performance 
of our method. 



C. Polytope packing 

By inspection of a packing of regular tetrahedra yielded 
by our numerical search during early phases of its devel- 
opment, we were able to construct a new transitive, peri- 
odic (p — 4) packing of tetrahedra with a higher density 
{(j) sa 0.8547) than previously reported [1]. This packing 
takes the form of a double lattice of bipyramidal dimers 
(the union of two face-sharing tetrahedra) . The packing 
has since been slightly improved to a closely related, but 
less symmetric packing with density w 0.8563 ITlj. 
In its current form, our search method is able to repro- 
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d 


-^highest 


(L) 

T 

highest 


(Alitor) 


(n) 


success rate 


2 


A2 


6 


27 


12 


100/100 


3 


D3 


12 


54 


40 


100/100 


4 


D4 


24 


132 


118 


98/100 


5 


D5 


40 


163 


331 


94/100 


6 


Ee 


72 


225 


928 


64/100 


7 


Et 


126 


597 


2729 


66/100 


8 


Eg 


240 


511 


6988 


55/100 


9 


Ag 


272 


350 


15604 63/100 


10 


Aio 


336 


438 


32203 28/100 


11 


All 


438 


549 


73766 10/100 



TABLE III: Results of PDC searches for lattice packing with 
high kissing number in dimensions d = 2, ... 11. For each 
dimension, 100 runs from random initial conditions were per- 



formed with a target coordination rta 



the high- 



est coordination number known for a lattice of that dimen- 
sion, Ahighost 'QM- The runs were limited to 5000 iterations, 
and the number of converged runs is quoted in the right-most 
column. The mean number of difference map iteration in 
converged runs was (Mter), and the mean number of relevant 
exclusion constraints used was (n). 
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FIG. 4: The course of a sample run searching for dense peri- 
odic packings {p — 4) of unit edge-length regular tetrahedra, 
showing ^Ai, a measure of the total interpenetration be- 
tween tetrahedra in the "concur" estimate (blue, defined in 
(|33|l '). and e^, the squared distance between the "divide" and 
"concur" estimates (purple, shifted up for clarity), both on a 
logarithmic scale. The density target for the search is started 
at ^target =0.75 aud adjusted when the search is converged 
on a solution (vertical red lines) to ^target ~ 0.82 (at itera- 
tion 15751) and then to (^target = 0.8563 (at iteration 15898). 
Each iteration took 14 millisecond on average on a single 3 
GHz CPU. 



primitive pentatopes 



where 



lattice 



where 



Ki = conv{ri,r2,r3,r4,r5} 
K2 = conv{r2,r3,r4,r5 re} 

Ki = t-K2 
n = 75(1,1,1,1) 
r2 = (3, -1,-1,-1) 
r3 = (-1,3, -1,-1) 
r4 = (-1,-1,3,-1) 
r5 = (-1,-1,-1,3) 
re = -^(1,1, 1,1) 
t= 1(-7,1,3,3)-^(1, 1,1,1 



A : 



Mo 



/-6 



10 

-8 -4 



I ^5 

^ 4 



-7 

/ 2 2 

2 2 
1 1 

3 3 



-7 
2 2 \ 

2 2 
1 1 

3 3 



-6 2 
4 8 
9 
9 



8 

-7 

-3/ 



TABLE IV: Coordinates of the densest pentatope packing 
discovered by the PDC search {(j) = 4 vol(A'i)/ det(Mo) = 
128/219 ^ 0.5845). 



takes the form of a double lattice of dimers (a dimer here 
is the union of two cell-sharing pentatopes). This struc- 
ture, composed of a repeating unit of two oppositely ori- 
ented dimers, repeatedly came up as the densest in de 
novo PDC searches with p = 4 and p = 8 pentatopes in 
the unit cell, whereas searches with intermediate values 
of p yielded sparser packings. We subsequently refined 
the packing with a restricted search where the dimer was 
taken as the basic particle. 

Note that the density reported is slightly lower than 
that of the densest known packing of four-dimensional 
spheres {(j) = 7r^/16 ~ 0.6169). It remains to be deter- 
mined whether this is the case because the optimal pack- 
ing density of pentatopes is smaller than that of spheres 
or because the dimer double lattice is suboptimal. The 
vertex coordinates of the four primitive pentatopes and 
the generating matrix of the lattice are given in Table 
lYl 



VI. CONCLUSION 



duce this densest known packing reliably (fifteen out of a 
hundred runs converged within the iteration limit), and 
Figure |3] shows the results of a sample run converging to 
this packing. 

For the problem of packing regular four-dimensional 
simplices (pentatopes) in four-dimensional Euclidean 
space, we report a new packing discovered by our search 
method (Figure [5|). This packing, with density (j) — 
128/219 w 0.5845, is, to our knowledge, denser than any 
previously reported packing of regular pentatopes. Like 
the densest known tetrahedron packing, this packing also 



In this article we report on the development of PDC, 
a novel, constraint-based method for discovering dense 
periodic packings through de novo numerical searches. 
We lay out the principles of the method and demon- 
strate its application for selected problems. In addition 
to the dense packing of regular tetrahedra reported in 
Ref. [l|, we also discover a new dense packing of regu- 
lar pentatopes using the PDC method. We also use the 
method to numerically recover the lattice sphere pack- 
ings of highest known density and highest known kissing 
number in a range of dimensions, providing empirical ev- 
idence of their optimality. 
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FIG. 5: The top figure shows a two-dimensional cut through 
the densest known packing of tetrahedra. The plane of the 
cut is parallel to the bases of the bipyramidal dimers. Trian- 
gular sections from dimers of one orientation (red) and from 
dimers of inverted orientation (blue) are visible. The bot- 
tom figure shows a three-dimensional cut through the densest 
known packing of pentatopes. The cut is taken parallel to the 
bases of the pentatope dimers, and tetrahedral sections from 
the two dimer orientations (red and blue, again) are visible. 



In developing the PDC scheme, we adapt the D — C 
framework to periodic systems. PDC retains the mind- 
set of the traditional D — C approach of Ref. [l^, but 
generalizes its formalism in a few ways. We introduce 
an expanded configuration space parameterized by lin- 



ear combinations of the original parameters, such that 
these new parameters over-determine the configuration. 
Therefore, by contrast with the traditional construction, 
where new parameters are, specifically, redundant copies 
of original parameters and concurrence is described by 
the equality of all copies of a given original parameter, 
here we allow concurrence to be described by a general 
linear relation. With this generalization, we can treat 
the periodic images of a particle as "replicas" of the par- 
ticle, even as they are related by a lattice vector instead 
of being identical. Thus, the variables describing the pe- 
riodic repetition of the configuration, namely the lattice 
vectors, are not imposed as constants or adjusted in ded- 
icated steps. Instead, due to the projection formulation 
of the dynamics, the unit cell variables that minimize the 
change to the configuration are determined at each iter- 
ation. These variables are treated on the same footing as 
particle positions and orientations and are optimized as 
aggressively. 

Additionally, we develop a displacement-minimizing 
overlap resolution algorithm for the convex hulls of two 
sets of points in M''. We use this algorithm to implement 
the projection to the exclusion constraint in the case of 
polytopal particles. 

Unlike Monte Carlo simulations, which explore the 
physical optimization landscape using stochastic moves, 
a PDC search uses a deterministic map in an expanded, 
non-physical configuration space. As such, it is useful 
when interest lies more in discovering optimal configu- 
rations and less in discovering the physical pathways to 
such configurations. However, introducing non-physical 
dynamics has been observed to be important in over- 
coming dynamical stagnation Q. The projection-based 
dynamics make PDC particularly well-suited in problems 
with hard constraints, such as hard particle packing, or 
with step potentials, which prohibit the use of gradient 
information. 

While no direct comparison has been made between 
the performance of PDC and Monte Carlo searches in 
the case of periodic packing problems, difference map and 
D — C methods in the case of other problems have been 
shown to perform better than or on a par with specialized 
and general-purpose methods [l8l - l20l .l29| . The generality 
of the PDC scheme and its demonstrated ability to dis- 
cover dense packings in a variety of settings indicate its 
utility as a general method for conducting de novo nu- 
merical searches and as a possibly attractive alternative 
to conventional methods [30|. 

Y. K. acknowledges N. Duane Loh for valuable dis- 
cussions. This work was supported by grant NSF-DMR- 
0426568. 
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